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Fortran 77 programs for the computation of modified Bessel functions of purely imaginary order 
are presented. The codes compute the functions Ki a (x), Li a (x) and their derivatives for real a 
and positive x; these functions are independent solutions of the differential equation x 2 w" + xw' + 
(a 2 — x 2 )w = 0. The code also computes exponentially scaled functions. The range of computation 
is (x, a) £ (0,1500] X [—1500,1500] when scaled functions are considered and it is larger than 
(0, 500] X [—400, 400] for standard IEEE double precision arithmetic. The relative accuracy is 
better than 10" 13 in the range (0, 200] X [-200, 200] and close to 10~ 12 in (0, 1500] X [-1500, 1500]. 

Categories and Subject Descriptors: G.4 [Mathematics of Computing]: Mathematical software 
General Terms: Algorithms 

Additional Key Words and Phrases: Bessel functions, numerical quadrature, asymptotic expan- 
sions 



1. INTRODUCTION 

These routines compute the functions K ia (x) and L ia (x), which constitute a nu- 
merically satisfactory pair of independent solutions of the modified Bessel equation 
of imaginary order: 

X 2 w" + xw' + (a 2 -x 2 )w = 0. (1) 

The routine DKIA computes Ki a (x) and its derivative. DLIA computes Li a (x) 
and its derivative. Both routines require that a method for computing Airy func- 
tions of real variable is available. The routines use Algorithm 819 [2] • 

The algorithm combines different methods of evaluation in different regions. 
Comparison between the different methods (see accompanying paper PP), together 
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with the use of the Wronskian relation, has been used to determine the accuracy 
of the algorithm (see Section [21 . 

The algorithm computes either the functions K ia (x), L ia (x) and their derivatives 
or scaled functions which can be used in larger regions than the unsealed functions. 
The scaled functions are defined in the accompanying paper pQ (Section 3, Eqs. 
(43) and (44)) and in the comment lines of the algorithms. The relative accuracy 
of the algorithms is (except close to the zeros of the functions) better than 1CP 13 
in the (x, a)-region Di = (0,200] x [-200,200], better than 5 x 10~ 13 in D 2 = 
(0,500] x [-500,500] and close to 10~ 12 in D 3 = (0,1500] x [-1500,1500]. Both 
scaled and unsealed functions can be computed in D± but only scaled functions can 
be evaluated in the whole domains D 2 and D% without causing overflows/underflows 
for typical IEEE machines. See the accompanying paper f]], Section 3. 

2. REGIONS OF APPLICATION OF THE METHODS AND ACCURACY. 

Because the scaled and unsealed functions are even functions of a, we can restrict 
the study to a > 0. 

In order to determine the region of applicability of each method, we have com- 
pared all the methods available with the non-oscillating integral representations 
in the accompanying paper [I], Section 2.4, which are expected to be valid in all 
the (x,a) plane except close to a = x. This is so because, as explained in [3], the 
integration paths become increasingly non-smooth as one approaches the transition 
line a = x (in this region, the uniform asymptotic expansions in the accompanying 
paper, Section 2.3, are the best option). The remaining approaches have little or no 
overlap in their regions of validity, except for the continued fraction method for the 
computation of Ki a (x) and K' ia (x) and the asymptotic expansions for large x; in 
fact, the validity region of the CF-method completely covers that of the asymptotic 
expansion; therefore, the CF method is preferred over the asymptotic expansion for 
large x for the computation of Ki a {x) and its derivative. 

As we will next describe, for not too large values of a and x there are at least 
two alternative methods of computation in any region of the (x, a) plane. This is 
important in order to determine the accuracy reachable by the algorithm. A second 
validation of the methods is provided by the Wronskian relation 

Ki a {x)L' ia {x) - K' ia (x)L ia (x) = 1/x, 

which is also satisfied by the scaled functions Ki a (x), K' ia (x), Li a (x), L' ia (x). 

Of course, when two different approaches are available within a given accuracy 
in a same region, we choose the fastest method of computation. By order of speed, 
from fastest to slowest we can list the different methods of computation described 
in £Q as follows: 

(1) Series expansions and continued fraction method. 

(2) Asymptotic expansions. 

(3) Integral representations. 

Figures 1,2 show the regions in the plane (x, a) where the different methods are 
used in the routines DKIA and DLIA, respectively. In the first comment lines of 
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the codes, the explicit equations for the curves separating the different regions are 
given. 
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Figure 1. Regions of applicability of each method for the Ki a (x), K' ia (x) functions. UAE: 
uniform Airy-type asymptotic expansion; I: integral representations; S: power series; CF: for 
continued fraction and AE: asymptotic expansions. 
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Figure 2. Regions of applicability of each method for the Li a (x), L' ia (x) functions. UAE: 
uniform Airy-type asymptotic expansion; I: integral representations; S: power series; CF: for 
continued fraction and AE: asymptotic expansions. 

Both the direct comparison with integral representations as well as the Wron- 
skian check show that the unsealed and scaled functions can be computed with an 
accuracy better than 10~ 13 in the region (x, a) e (0,200] x [0,200] without using 
integral representations. This accuracy was found to be in a good compromise with 
efficiency We expect that function values outside this range will not be needed very 
often. As the range considered is increased, the accuracy decreases mildly, being 
better than 5 x 10~ 13 in (0,500] x [0,500] and close to 10~ 12 in the full range of 
computation (0, 1500] x [0, 1500]. Of course, near the zeros of the functions (there 
are infinitely many of them in the oscillatory region x < a) relative accuracy loses 
meaning and only absolute accuracy makes sense. 

3. TIMING 

Figure 3 shows the CPU time spent by the code in the regions (x, a) G (0, L] x [0, L]; 
these CPU times refer to a Pentium-II 350-Mhz processor running under Debian- 



4 • Amparo Gil, Javier Segura and Nico M. Temme 

Linux; the compiler used was the GNU Fortran compiler g77. 
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Figure 3. CPU times in (0, L] X [0, L] as a function of L. 

It is apparent from the figure that different methods with different speeds appear 
as L grows. For small L, only series (and the CF method for Ki a {x)) are needed, 
which are the fastest methods. As soon as asymptotic expansions come into play 
the CPU times tend to increase more rapidly. For L larger than 200 the effect of 
the use of quadratures becomes prominent. 

Tables 1 and 2 show the average time spent by each method of computation in 
microseconds, both for DKIA and DLIA routines. Also, the percentage of use of 
each method is shown. These data are shown for two different regions of computa- 
tion: the "10~ 13 accuracy region" ((0,200] x [0,200]) (Table 1) and the full range 
permitted by the routines ((0, 1500] x [0, 1500]) (Table 2). 
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Table 1 . Percentage of use and average CPU times of each method in the region (0, 200] X [0, 200] . 
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Table 2. Percentage of use and average CPU times of each method in the region (0, 1500] X 
[0, 1500]. 

The algorithm has been tested in several computers and operating systems (Pen- 
tium II PC under Debian Linux, Pentium IV laptop under Windows XP, Pentium 
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IV PC under Windows XP and Red Hat Linux and Solaris 8 workstation) and 
several compilers (g77 and f77 for Linux and Windows, Digital Fortan and Lahey 
Fortran for Windows and the Sun Fortran 95). 



4. COMPARISON WITH EXISTING SOFTWARE 

In spite of the many applications of the modified Bessel functions of imaginary 
order, the only previously published algorithm, as far as the authors know, is the 
code by Thompson and Barnett |3] (TB algorithm from now on). This code is 
intended for the computation of Coulomb wave functions of complex order and 
arguments, which have as a subset modified Bessel functions of complex orders. It 
is not surprising that, as we later discuss, the TB algorithm has limitations when 
computing the Ki a {x) and Li a {x) functions: this is a more general algorithm than 
ours and it was specially designed for real parameters, although it can be applied 
for complex parameters too. Apart from this code, there are no available algorithms 
in the public domain, as can be judged by browsing at GAMS (Guide to Available 
Mathematical Software: http://gams.nist.org), which are able to compute Ki a (x), 
Li a (x) and their derivatives. 

As Figure 4 shows, the TB algorithm (program COULCC available from CPC 
library) is inaccurate for imaginary orders and suffers from overflow and convergence 
problems. We have tested the program COULCC 0] against our algorithms. Figure 
4 shows the comparison. At the lighter points the flag for detecting errors in 
COULCC was different from zero, showing that the algorithm failed to converge 
or suffered overflows. The rest of points show discrepancies with our code, which 
persist for 10 -8 accuracy and lower. Of course, at the lighter points all the accuracy 
is lost and the problems persist for larger values of x and a. 
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Figure 4. Regions where the TB algorithm fails due to poor convergence or overflows (light 
shaded points) and points where accuracy is lost (dark shaded points). The discrepancy is greater 
than 10 — 8 in the left figure and greater than 10 — 13 in the right figure. 

Therefore, our algorithm seems to be the first accurate and efficient code which is 
capable of computing modified Bessel functions of imaginary order with an accuracy 
close to full double precision in a wide domain. 
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